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ABSTRACT 

o' 
o , 

The observed abundance of high-redshift galaxies and clusters contains precious 
information about the properties of the initial perturbations. We present a method to 
compute analytically the number density of objects as a function of mass and redshift 
■ for a range of physically motivated non-Gaussian models. In these models the non- 

l/-) ■ Gaussianity can be dialed from zero and is assumed to be small. We compute the 

probability density function for the smoothed dark matter density field and we extend 
{Sj . the Press and Schechter approach to mildly non-Gaussian density fields. The abundance 

of high-redshift objects can be directly related to the non-Gaussianity parameter and 
\Q \ thus to the physical processes that generated deviations from the Gaussian behaviour. 

Even a skewness parameter of order 0.1 implies a dramatic change in the predicted 
abundance of z > 1 objects. Observations from NGST and X-ray satellites (XMM) can 
be used to accurately measure the amount of non-Gaussianity in the primordial density 
field. 

Q-f 

q . Subject headings: Cosmology: theory — galaxies: clusters 



1. Introduction 



In the standard inflationary model for triggering structure formation in the Universe, there 
are precise predictions for the properties of the initial fluctuations: they are adiabatic and follow 
a nearly-Gaussian distribution, with deviations from Gaussianity which are calculable, small and 
generally dependent on the specific inflationary model (Falk et al. 1993; Gangui et al. 1994; Gangui 
1994; Wang & Kamionkowski 1999; Gangui & Martin 1999). Because of the smallness of such 
deviations from the Gaussian behaviour, in most theoretical predictions one simply assumes that the 
primordial density field has exactly random phases. Consequently, intrinsic temperature fluctuations 
in the Cosmic Microwave Background (CMB) are commonly treated as being Gaussian and the same 
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assumption is made in most analyses of Large-Scale Structure (LSS) of the Universe. Besides this 
'standard' nearly-Gaussian model for generating cosmological structures, based on the amplification 
of quantum fluctuations of the same scalar field which drives the inflationary dynamics, there exist 
alternative models for the origin of fluctuations which predict stronger deviations from the random- 
phase paradigm. Still within the context of inflation multiple scalar field models can give rise to 
non-Gaussian perturbations of either isocurvature or adiabatic type (Allen et al. 1987; Kofman 
& Pogosyan 1988; Salopek et al. 1989; Linde & Mukhanov 1997; Peebles 1999a,b; Salopek 1999). 
Alternatively, cosmological defect scenarios (e.g. Vilenkin 1985; Vachaspati 1986; Hill et al. 1989; 
Turok 1989; Albrecht & Stebbins 1992) generally predict non-Gaussian initial conditions. 

The observed abundance of high-redshift cosmic structures contains important information 
about the properties of initial conditions on galaxy and clusters scales. CMB observations will put 
constraints on the nature of the initial conditions (e.g. Pen & Spergel 1995; Hu et al. 1997; Verde 
et al. 1999; Wang & Kamionkowski 1999; Wang et al. 1999), but will be severely affected by the 
presence of noise and foregrounds (e.g. Knox 1999; Tegmark 1998; Bouchet k, Gispart 1999) and 
probe rather large scales. 

The Gaussian assumption plays a central role in analytical predictions for the abundance 
and statistical properties of the first objects to collapse in the Universe. In this context, the 
formalism proposed by Press & Schechter (Press & Schechter 1974), with its later extensions and 
improvements (Peacock &: Heavens 1990; Bond et al. 1991; Cole 1991; Lacey &: Cole 1993) has 
become the 'standard lore' for predicting the number of collapsed dark matter halos as a function 
of redshift. However, even a small deviation from Gaussianity would have a deep impact on those 
statistics which probe the tails of the distribution. This is indeed the case for the abundance of 
high-redshift objects like galaxies at z > 5 or clusters at z ~ 1 which correspond to high peaks, 
i.e. rare events, in the underlying dark matter density field. Therefore, even small deviations from 
Gaussianity might be potentially detectable by looking at the statistics of high-redshift systems. 

The importance of using the mass-function as a tool to distinguish among different non- 
Gaussian statistics for the primordial density field, was first recognized by Lucchin & Matarrese 
(1988), Colafrancesco et al. (1989) and, more recently, by Chiu et al. (1998), followed by Robinson 
& Baker (1999), Robinson et al. (1999a,b), Koyama et al. (1999), Willick (1999), Avelino & Viana 
(1999). To make predictions on the number counts of high-redshift structures in the context of 
non-Gaussian initial conditions, a generalized version of the PS approach had to be introduced. 

Such a generalization of the PS formalism has been tested successfully against N-body simu- 
lations (Robinson & Baker 1999), but, from a theoretical point of view, it suffers from the same 
problems of the original Press & Schechter (PS) formulation: it cannot properly account for the 
so called cloud-in- cloud problem, i.e. the constraint that bound systems of given mass should not 
be incorporated in larger mass condensations of the same catalog 4 . Even more important, most of 
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the PS generalizations to non-Gaussian models proposed in the literature do not properly take into 
account the dependence of the fluctuation field on the smoothing scale. 

In this paper we obtain an analytic prediction for the number of dark matter halos as a function 
of redshift, within the hierarchical structure formation paradigm, following the PS formalism. The 
strength of our method is that we are able to properly take into account the smoothing-scale, 
or mass, dependence of the probability distribution function of the primordial density contrast. 
Obtaining analytical results in this context is extremely important. Direct simulations of non- 
Gaussian fields are generally plagued by the difficulty of properly accounting for the non-linear 
way in which resolution and finite box-size effects, present in any realization of the underlying 
Gaussian process, propagate into the statistical properties of the non-Gaussian field. Moreover, 
finite volume realizations of non-Gaussian fields might fail in producing fair samples of the assumed 
statistical distribution, i.e. ensemble and (finite-volume) spatial distributions might sensibly differ 
(e.g. Zel'dovich et al. 1987). This problem, of course, becomes exacerbated and hard to keep under 
control in so far as the tails of the distribution are concerned. Thus, in looking for the likelihood 
of rare events for a non-Gaussian density field, either exact or approximate analytical estimates 
should be considered as the primary tool. 

When compared with N-body simulation results the striking feature of the standard PS algo- 
rithm is that it works extremely well in predicting properties of highly non- linear objects such as 
the DM haloes, even being based on linear theory and relatively simple assumptions. This agree- 
ment has been proven for Gaussian initial conditions, and one might wonder whether such a feature 
also extends to the non-Gaussian case, where additional non-linear couplings arise and non-linear 
gravitational corrections can be important at larger scales, unlike the Gaussian case. This exten- 
sion of the PS prediction will be tested against N-body simulations in a subsequent paper. For the 
present application, however, the primordial non-Gaussianity is assumed to be small and we can 
safely assume that departures from the PS prediction are negligible. 

The outline of the paper is as follows: in Section 2 we propose a parameterization of primordial 
non-Gaussianity that covers a wide range of physically motivated models whose non-Gaussianity 
can be dialed from zero (the Gaussian limit). We will then relate analytically the non-Gaussianity 
parameter to the number of high-redshift objects. This step involves the generalization of the 
PS formalism to non-Gaussian initial conditions, which in turn requires an expression for the 
probability density function of the non-Gaussian mass-density field as a function of the filtering 
radius. In Section 3, we calculate the probability density function for the density fluctuation field, 
smoothed on the mass scale typical of high-redshift objects. Since the non-Gaussianity is expected 
to be small, we can calculate the probability density function by expanding the 'cumulant generator' 
in powers of the non-Gaussianity parameter and keeping only linear order terms. In this context 
we are able to obtain an analytic expression for the probability density function of the smoothed 
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density field. We generalize the PS approach to non-Gaussian density fields in Section 4, where we 
give an analytical expression for the comoving mass- function of halos formed at a given redshift, 
also accounting for the cloud-in-cloud problem. Section 5 discusses the association between halos 
and high-redshift galaxies. The PS formalism predicts the number of dark matter halos, but the 
correspondence to the observed number of galaxies is not necessarily simple. Finally in Section 6 
we summarize our main results and draw our conclusions. 

2. Parameterization of primordial non-Gaussianity 

Because of the infinite range of possible non-Gaussian models, we consider here two models 
with primordial non-Gaussianity, whose amplitude can be dialed from zero (the Gaussian limit). 
Following Verde et al. (1999) we consider models in which either the density contrast (model A) 
or the gravitational potential (model B) contains a part which is the square of a Gaussian random 
field. The physical motivation for this choice is that such a non-Gaussianity may arise in slow-roll 
and/or nonstandard (e.g. two-field) inflation models (Luo 1994; Falk et al. 1993; Gangui et al. 
1994; Fan h Bardeen 1992); moreover both models may be considered as a Taylor expansion of 
more general non-Gaussian fields (e.g. Coles & Barrow 1987; Verde et al. 1999), and are thus fairly 
generic forms of mild non-Gaussianity. 

We note here that, since in the PS framework the evolution of perturbations is considered to be 
linear, this parameterization of non-Gaussianity is effectively equivalent to non-linear biasing acting 
on a truly Gaussian underlying field. The effect of biasing is to alter the clustering properties of the 
galaxies with respect to the dark matter by relating galaxy formation efficiency to the environment. 
We do not intend to investigate the clustering properties of high-redshift objects, we will instead 
focus on the probability density function (PDF) of the linear density contrast as a tool to predict 
the abundance of dark matter halos as a function of mass and formation redshift. 

2.1. Linear plus quadratic model for the density or the gravitational potential 

As a first step, we want to obtain the PDF for a model where either the primordial over-density 
field 5 = Sp/p, or the initial peculiar gravitational potential $ is represented by a zero-mean random 
field ip, given by the following local transformation T on an underlying Gaussian field <fi: 

V(x) = F[<j>] = a#x) + e(^ 2 (x) - <</> 2 » , (1) 

with a and e free parameters of the model. The homogeneous and isotropic Gaussian process <f> is 
assumed to have zero mean, ((f)) = 0, and power-spectrum to be specified later. Note that in 
the limit a — ► 0, ip is chi-squared distributed, while for e — ► one recovers the Gaussian case. As 
pointed out before, the model of eq. (1) can be thought as representing the first two terms of the 
Taylor expansion of more general non-Gaussian fields around the Gaussian limit. Cubic or higher 
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order powers of cf> in such an expansion should be thought as being of order e 2 or higher. We will 
come back to this point in Section 3.3. Since the overall normalization of ip is fixed observationally 
(e.g. by measurements of the power-spectrum), effectively there is only one free parameter in the 
model. In what follows we will take either a = 1 and e as the non-Gaussianity parameter or a = 
and e = 1 to obtain the \ 2 model. 

If there was no filtering to take into account or if the transformation (1) were true for a 
smoothed field, then the PDF for ip could be simply obtained as follows 

P(iP) = (S D (a<P + e(<p 2 - (<p 2 )) = ] d<pP(<P)8 D (acP + e{<p 2 - {<p 2 )) - ip) (2) 

where P{4>) is a Gaussian PDF with zero mean and variance (cp 2 ) and 5 D is the Dirac delta function. 
The above equation can be also thought as resulting from applying the Chapman-Kolmogorov 
equation: P(ip) = J d(pW(ip\(p)P{(p) where W(tp\cp) = 5 D {ip — T[<p\) is the transition probability 
from <p to ip (e.g. Taylor h Watts 2000). The integral over <p can be performed by writing: 

P(V) = J d<pP(<P)5 D {<P - jr-i^^L. , (3) 

which gives: 

P(iP) = {27T<0 2 > [a 2 + 46 {e{cp 2 ) + */>)] [S-W ~ £+&)] , (4) 

with 

£±(V0 = exp | -^^y [a 2 + 2e (e(0 2 ) +i/>)± a^a 2 + 4e(e((f> 2 ) + V)] j • (5) 

However, in order to apply the PS formalism, we need to know the PDF for the linear mass- 
density field as a function of the smoothing radius R and redshift z. To take into account the 
different evolution of modes inside the horizon, we assume that the transfer function acts on our 
primordial non-Gaussian field ^(x), as a convolution. We can easily account for the smoothing 
operation, the effect of the transfer function and the linear growth of perturbations, by writing 
(e.g. Moscardini et al. 1991): 

<fo(x, z) = D(z)5 R (x) = D{z) j d 3 yF R (\x - y|)^(y) , (6) 

where D(z) is the growing mode of linear perturbations, normalized to unity at z = 0, so that 
<5r(x) represents the mass density fluctuation field linearly extrapolated to the present time. The 
isotropic function Fr(|x|) can be specified by its Fourier transform, 

F R (k) = j d 3 y e-^F(|y|) = W R (k)T(k)g(k) , (7) 

where W R (k) is the Fourier transform of the assumed low-pass filter (e.g. a spherical top-hat filter), 
T(k) is the transfer function (normalized to unity for k — > 0), that we will later assume to take the 
adiabatic CDM form as in (Bardeen et al. 1986), with the modifications of Sugiyama (1995). 

The function g(k) completes the specification of our model: 
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• model A: g(k) = 1, if we assume that our non-Gaussian model applies directly to the 
primordial density field. 

• model B: g(k) = —^(k/Ho) 2 ^}^, if the same assumption is made for the gravitational 
potential (its precise form comes from solving the cosmic Poisson equation); here $7o m denotes 
the present day density parameter of non-relativistic - baryonic plus dark - matter. 

As widely discussed by Moscardini et al. (1991), the overall sign of the Gaussian to non-Gaussian 
mapping inherent in these models is a crucial parameter, which determines most of the non-linear 
dynamics. In our case, we are still free to choose the sign, through the actual choice of the free 
parameter e. We finally get 

6 R (yi) = a0 K (x) + e J d 3 yF R (\ x - y |)0 2 (y) - C , (8) 

where <^r(x) is just the smoothed underlying Gaussian field i.e. the convolution of <p with F, and 

C^e{ct> 2 ) J d 3 yF R (\ x-y |) (9) 

ensures that Sr has zero expectation value. 



3. PDF of the smoothed density fluctuations 

In this section we derive an approximate analytic expression, valid for small values of the non- 
Gaussian parameter e, for the PDF as a function of the smoothing radius (i.e. mass scale) and of 
the redshift of collapse. This is the key ingredient to develop the extension of the PS formalism 
to non-Gaussian fields. The reader mainly interested in the application of the PS approach to 
non-Gaussian density fields may omit reading this section at a first sitting. 

Once the filtering process has been taken into account, the density fluctuation PDF can be 
obtained through a functional or path-integral 5 (e.g. Ramond 1989) 

P(5r) = (S D (a<f> R (x) + e J d 3 yF R (\ x - y |)0 2 (y) - C - 5 R (x^J > (10) 

= J [Vc(>]V[cf>] j ^ exp iX (a<Mx) + * / d 3 yF R (\ x-y |)0 2 (y) - C - Sr(^ , 

where, in the second line, we used the integral representation of the Dirac delta function. The 
functional integration is over 4> configurations in real space weighted by the Gaussian probability 



5 The path-integral approach has been widely applied in the cosmological context and in particular to large-scale 
structure studies by e.g. Politzer & Wise (1984); Grinstein & Wise (1986); Matarrese, Lucchin & Bonometto (1986); 
Bertschinger (1987). 
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density functional 

P[<f>] =exp|-^ J d 3 y j d 3 #(y)/C(y,z)0(z)| / J exp { - ^ J d 3 y j d 3 ^(y)/C(y, 
which has been consistently normalized to unit total probability, J[T><j)]'P[(/)\ = 1. 



WW 
(11) 



Although the previous expressions for the PDF actually gives its form at redshift z = 0, one 
should keep in mind that the quantity P(5 R )d5 R is redshift-independent, as long as linear evolution 
applies. 

From now on we will use the following compact notation: 

J d 3 y j d 3 #(y)/C(y,z)0(z) = (12) 

The symmetric kernel K, is defined as the functional inverse of the two-point correlation function 
£tf>(y> w) = ^(|y - w|) of the field <j>, namely, 

j d 3 y/C(z,y)^(y,w) =5 D (z- w) . (13) 

By exploiting the fact that 

j d 3 yF R {\ x - y \)<p 2 (y) = j d 3 y j d 3 z<P(y)F R (\ x -y\)5 D (y- z)0(z) = (<f>, F R 5 D , <fi), (14) 

equation (11) becomes: 

p(x v _ / fe-^-^ exp {-1(0, /C, 0) + »Ae(& F R 6 D , 0) + ^ 0)} 

where we have defined the source functional 

J\{y) = A<xFr(| x — y |) (16) 

and introduced the notation 

(Jx,(f>) = J d 3 yj x (y)4>(y) . (17) 

The above functional integration can be performed analytically, by applying the so-called path- 
integral technique for composite operators (Cornwall et al. 1974; Hawking & Moss 1983). Let us 
briefly sketch the main steps of the procedure. 

We start by defining a new kernel: 

/C' A (z, y) = /C(z, y) - i2XeF R (\ y - x \)5 D (z - y) , (18) 
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which allows to write 



P{,) ti I ^ e -^«-^ c [[v^e-^'x^+KJ^) . (19) 



We then make the change of variable (under which the path-integral is left unchanged) 

<f>(z) - 0(z) - ^ d 3 y [/CI]" 1 (z,y)J A (y) (20) 

where [/C^J -1 satisfies the integral equation 

j A/C , (z,y)[/Cl]- 1 (y,w)=^(z-w). (21) 

We then easily get: 

P{5r ) = J ^ e -i*S R -iXC-l(j x ,[K' x ]-\j x ) { Det [^]/Det[X;]}- 1/2 , (22) 

having used the standard notation for functional determinants: 

{Det[K' x ]/Det[K\}-^ = im^i-^xM = exp (_I T rln [l - i2Xe(F R S D , XT 1 )] ) . 

(23) 

The reader unfamiliar with the functional notation can understand the last result as a generalization 
of the identity ln[detM] = Tr[lnM], which applies to any symmetric matrix M. Here 1 denotes 
the functional unit matrix, i.e. the Dirac delta function; the logarithm of a functional is defined as 
its series expansion and 

{F R 5 D ,K,- 1 )^ J d 3 wF R (\w- X \)5 D (y-w)IC-\w-z) = F R (\y- X \)^(y-z) . (24) 

Finally, the trace (Tr) of a functional G(y, z) is defined as / d 3 y J d 3 z5 D (y — z)G(y,z). 
Therefore, equation (23) becomes: 

{Det[K,' x }/Det[1C]}- 1 / 2 = exp |~ J d 3 y J d 3 z5 D (y - z) In [5 D (y - z) - i2XeF R (\y - x|)^(y - z)] 

(25) 

The PDF for the linearly evolved, smoothed density field takes the exact form: 

P(S R )dS R = J ^ e -^+W(A) d§R (26) 

where W(A) is called the cumulant generator, as its series expansion around A = defines the 
cumulants (or irreducible moments) of S R (x). Its exact form is: 



W (A) = -i\C-±jd 3 yjd 3 z! [ \WF R (\y\)[lC x ]- 1 (y,z)F R (\z\) 
+ < 5 D (y-z)ln[,5 D (y-z)- i 2A6F R (|y|)^(y-z)]| . 



(27) 
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where the A subscript indicates where the A dependence is hidden. Note that the x dependence 
has been eliminated by a mere translation of the origin. Equation (27) contains also the exact 
form for the generating function for the particular case where the original non-Gaussian field ip is 
chi-squared distributed: this is simply obtained by setting a = in the previous expression. Note 
also that the first-order term in the expansion of the logarithm precisely cancels the — iXC term, 
so that the condition (Sr) = is identically satisfied. 



3.1. Kernel inversion 

In order to solve the remaining integrals we need to find an expression for the functional 
[/C^] _1 . Let us start by finding an expression for the kernel /C, given that its inverse is just the 
auto-correlation function of the Gaussian field 4>, namely, 

j d 3 y/C(w - y)^(| y-z |) = 5 D (w - z) . (28) 

Fourier transforming both sides of eq. (28) we easily find (Politzer & Wise 1984; Bertschinger 
1987) 

K[k) = l/P^k) (29) 

where P<f,(k) is the power-spectrum of the underlying Gaussian field <j>. Therefore 
/C(w, y) = /C(| w — y |) = /C(r) reads 



If 11 f°° 

/C(r) = I ^ exp( * k • r) P4k) = ^L dkk2jo{kr) 



P<t>(k) 



(30) 



To find an expression for [/C^J 1 one would like to proceed in an analogous way. However, 
owing to the absence of translational invariance of /C', going to momentum space does not help. 
The technique we are then going to use is to expand [ZC^J -1 in powers of the non-Gaussianity 
parameter e, as follows: 

oo 

[/C^]- 1 ( y ,z) = ^(2^Ar^)(y,z) (31) 

n=0 

where the coefficients IZ^ are obtained recursively from 

7^)(y,z) = |dV^|y-w|)F i? (|w|)^- 1 )(w,z) (n > 0) , (32) 
with 7£(°)(y, z) = ^(|y — z|). Similarly, in Fourier space we have 

7^)(k,k') = P (t> {k) j ^0i? i? (^(«- 1 )( q + k,k') , (33) 
with £(°)(k,k') = ^^(^(k + k'). 
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From equation (33) we find: 

K {1 \k,k') = P (f) (k)P lP {k')F R (\k + k'\) (34) 

and, for n > 2, 

£W(k,k') = P ^k)P^(k') J f^^F R { qi )...F R {q n ^) (35) 

X F R (\qi + ■■■ + q n _i + k + k'|)P^(|qi + k|) • • • P^(| qi + • • • + q n _i + k|) . 

Note that, while in principle the expressions for K! , and should be symmetrized, for the 
purpose of calculating the PDF or the cumulants, this operation is not needed. 



3.2. Cumulant generator 

At this point we are ready to provide an expression for the cumulant generator, by first 
expanding it in powers of A. We have: 

n=2 

where jJL n , R denotes the cumulant of order n of the smoothed density contrast S R . 

The variance, the skewness and the kurtosis of the smoothed non-Gaussian density field are, 
respectively 6 , 

a 2 f°° ~ 
= o\ = (5 2 R ) = ^2 J o dkk 2 FUk)P4,(k) 



+ 



2 roo poo pi ____ I 

J dk^P^k) J dk'k' P^k') J dn F 2 R {^k 2 + k> 2 + 2kk'n) , (37) 



o 2 roo ^ poo ____ pi ____ i 

= <A) = -^J dktfFnWPtik) J dk'k' F R (k')P4k') J^d f iF R (^k 2 + k' 2 + 2kk> f i) (38) 

/d 3 ki f rfik-> f rfiki ~ ~ ~ 

Wfj (2^7 7 ^P (fci)P0(fc2)^(^(|k 1 -k 2 |)^(|k 2 -k 3 |)F ii (|k3-k 1 |). 

6 The integral over k in the sub-leading term of the variance can diverge for certain choices of the power-spectrum of 
the underlying Gaussian field and transfer function. In these cases one should bear in mind that any physical process 
originating the underlying field will necessarily provide the <j> power-spectrum with both infrared and ultraviolet cutoffs 
(the present-day horizon and the reheating scale, respectively, in the case of inflation-generated perturbations). In 
case this contribution dominates over the leading-order contribution even for e < 1 one can postulate the existence 
of a term proportional to e 2 in eq. (1) that will cancel out the sub-leading term in the variance. Alternatively eq. 
(1) and eq. (37) can be renormalized by choosing a and A$ such that ji2,R — fJ^ R when e«l. Our calculations still 
apply in this case provided one interprets e as e/a 2 . 
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and 



H,, R (5 R ) - 3{5lf 

/d^k-i f d 3 kn f d^ko ~ ~ ~ 

(2^/(2^/ j^P^P^P^ + k^Fnik^F^Fnih) (39) 

x F R (|k 1+ k 2 + k 3 |) + 96 e 4 / ^P^WW+WP+M 
x Ffldki - k 2 |)F fl (|k 2 - k 3 |)F R (|k 3 - k 4 |)Fi ? (|k 4 - ki|) . 

More in general, the cumulants are made of two contributions: a leading term of order e n ~ 2 
plus a subleading term of order e", 

where 

»% = 2^n\ a 2 J ^ / ^F R ( kl )Mk 2 )n^\^ (41) 

and 

M?i = 2^ 1 (n-l)!| ^ • • •/ • • • Wn)^(|ki-k 2 |) • • • Ffldk^x-I^D^dlc-kxl) 

(42) 

Note that, as anticipated, the sign of the parameter e, for a given model, fully determines the 
sign of the skewness as well as of all the other odd-order cumulants. 

The cumulants for the filtered chi-squared model are immediately recovered by taking a = 

(2) 

and e = 1 in the previous expressions, so that fin,R = ^ n R - It is worth noticing that eq. (42) 
supplies all the cumulants of a chi-squared field smoothed on scale R. This is an interesting result 
for models such as that recently proposed by Peebles (1999a,b), where non-Gaussian isocurvature 
fluctuations are obtained with a chi-squared distributed density field. 



3.2.1. Linear order in e 
To first order, the skewness depends linearly on e, therefore let us define 

Ss,R = eSi% = e$ R /{$ R )\ (43) 

where 

2 roc 

$R = ^2 I dkk 2 F 2 R (k)P^k) (44) 
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and 

/4!i = (4) = ^4 y dk^Fni^P^k) J dk'k' Fr^P^X) J dfi Fr(\J fc 2 + k' 2 + 2kk'n) 

(45) 

To this order in e there is no contribution from higher order moments. On very large scales, 
where the transfer function becomes unity, and for a power-law <p power-spectrum (P^ = A^k n or 
P^ = A^k n ~ 3 for model A or B, respectively), this leading-order skewness parameter S^j? becomes 
scale-independent, thus mimicking the behaviour induced by the gravitational instability (Peebles 
1980). For model A the variance, the skewness fJ^ R an d skewness parameter Ss t R are plotted vs. 
the radius R and the mass M in solar masses, M oc R 3 , in Figure 1 with the following assumptions: 
a = 1, top-hat filter in real space, transfer function as in (Sugiyama 1995), with baryonic matter 
density parameter O 6 = 0.015/i~ 2 , h = 0.65, f^om = 0.3, and scale-invariant primordial power- 
spectrum P^ oc k, normalized so that the present day r.m.s. fluctuation on a sphere of 8 hr 1 
Mpc is as = 0.99 (Viana Sz Liddle 1998), which simultaneously allows to best fit the local cluster 
abundance and the COBE data (e.g. Tegmark 1996). In what follows we will always assume that 
the approximation a\ = applies on all considered scales. As detailed below, at z = 0, for 
e < 0.01 in model A we will recover the standard PS mass-function on clusters scales, so that the 
relation between the observed abundance of clusters and the value of a§ keeps unchanged, in spite 
of our non-Gaussian assumption. 

To deal with model B we take P^ oc k~ 3 and, as above, we normalize the mass-density variance 
as in model A, with the same choice of cosmological parameters and transfer function. The variance, 
the skewness and skewness parameter S^ R are plotted vs. the radius R and the mass M in 
solar masses, M oc R 3 , in Figure 2 with the same assumptions as for model A in Figure 1. The 
variance for model B to linear order in e is identical to the variance for model A, but the skewness 
and skewness parameter are different: in particular notice that S^jj for model B has the opposite 
sign of e and, for a given value of e, its amplitude is many orders of magnitude smaller than for 
model A. It is important to realize, that both signs for e are generally allowed. In the inflationary 
case, both the sign and the magnitude of e are related to the inflationary slow-roll parameters ei n fl 
and r/infl (Gangui et al. 1994; Gangui 1994; Wang &: Kamionkowski 1999; Gangui & Martin 1999). 
In this model, the PS mass-function at z = is practically recovered on all scales, for all |e| < 200. 



3.2.2. Quadratic order in e 

For small deviations from Gaussianity, a first-order expansion in e is a valid approximation. 
We give here the expressions for the relevant quantities also to second order, but in any application 
we will neglect second or higher-order corrections. 

When expanding the cumulant generator to second order in e we obtain that the variance has 
a contribution oc e 2 as in eq. (37), the skewness remains the same as to first-order in e, the next to 
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leading term being oc e 3 , but there is a non- vanishing contribution to the kurtosis oc e 2 , namely 

x F i ?(fc 1 )F i? (A: 2 )F R (A:3)F R (|k 1 + k 2 + k3|) (46) 

For both models, the leading-order kurtosis parameter S^r = e 2 /j,^ r / (fi^ R ) 3 becomes scale- 
independent if the power-spectrum is a power-law; therefore on large scales, where the transfer 
function is unity, S^r becomes scale-independent. 

3.3. Higher order non-Gaussian contributions 

What happens to our results if we allow for higher order terms in the original definition of 
our primordial non-Gaussian field ipl Let us modify, for instance, our definition by adding a cubic 
term, as follows, 

V(x) = «0(x) + e(</> 2 (x) - (<f> 2 )) + eW(x) , (47) 

where f3 is a new independent parameter, which we assume to be of order unity. Such a cubic 
term might be easily accounted for in the functional integral approach, by using the so-called 
'integration by parts' relation (e.g. Ramond 1989). For the purpose of the present paper, however, 
we can account for such a cubic term through the modifications it induces in the lowest-order 
cumulants. The sub- leading term of the variance would in fact be modified by the addition of 

= )A% (48) 

which therefore appears as a 'renormalization' of the leading-order variance. The sub-leading 
skewness would also get an extra contribution of the same order e 3 , which we will not write here. 
The kurtosis, instead, would be modified already to leading order, gaining the extra piece 

<?A$ R = \<?k^% , (49) 

which still has the nice feature of appearing as a renormalization of the previously calculated 
leading-order kurtosis. 

At this point we can derive a general 7 and self-consistent approximation to the cumulant 
generator (and hence for both the differential and cumulative probability) up to order e 2 : 



mA)=-^a-4^ +e2 



(50) 



7 The key assumption that allowed us to obtain such a general treatment of a mildly non-Gaussian field is that it 
can be expanded as a local functional of an underlying Gaussian field. On the other hand, any form of non-locality 
such that it can be expressed as a convolution in real space can also be handled in a similar way, by suitably modifying 
our definition of the function Fr(x). 
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This is the explicit expression for W(A); no further terms could appear from higher-order non- 
Gaussianity to this order in e. Substituting this expansion in eq. (26) one gets an approximate 
expression for the PDF as an integral over A, which is valid for suitably small values of the non- 
Gaussian parameter e. 

For the rest of our calculations, however, we will assume that departures from Gaussianity are 
small and therefore we will retain only linear-order terms. 



To obtain the abundance of dark matter halos as a function of filtering radius R (or mass 
M oc i? 3 ) and redshift of collapse z c , one should first obtain the conditional probability that the 
density contrast equals the threshold for collapse <5 C , when filtered on scale R, provided it is below 
it on any larger scale R' . In the Gaussian case, and for sharp-k-space filter, this problem has 
been solved by a number of authors (Peacock & Heavens 1990; Cole 1991; Bond et al. 1991) by 
rephrasing it in terms of the problem of barrier first-crossing by a Markovian random walk. In 
the non-Gaussian case and/or for other types of filters, different techniques might also be useful. 
In particular, an alternative formulation, originally proposed by Jedamzik (1995) and successively 
implemented by Yano, Nagashima & Gouda (1996), Nagashima Sz Gouda (1997) and Lee & 
Shandarin (1998), allows to reduce the problem to the solution of the integral equation, 



where p~o m is the present-day mean density of non-relativistic matter, P(> 5 c \z c ,M) is the 
probability that 5m lies above the threshold 5 C (i.e. the fraction of volume where this happens) at 
a given redshift z c and n(M, z c ) is the required comoving mass-function for halos of mass between 
M and M + dM which formed at z c . The function P{M\M') denotes the conditional probability 
of finding a region with mass M overdense by 8 C or more, given that it is included in an isolated 
region of mass M' (> M). 

Let us start by first finding an expression for the L.H.S. of this equation, i.e. for the level- 
crossing probability P(> S c \z c ,M). 

In looking at this problem, it is convenient to think of the density fluctuation as being time- 
independent while giving a redshift dependence to the collapse threshold S c (z c ) = A c (z c )/D(z c \ 
^0m,^0A); A c is the linear extrapolation of the over-density for spherical collapse: it is 1.686 in 
the Einstein-de Sitter case, while it slightly depends on redshift (see Figure 3) for more general 
cosmologies (e.g. Kitayama & Suto 1996); Qqa denotes the closure density of vacuum energy today. 
Here D(z \ £l 0m ,Q x) denotes the general expression for the linear growth factor, which depends 
on the background cosmology: the redshift dependence of D(z c | f2om,^oA) _1 is shown in Figure 3 
for three different cosmological models. 



4. Press-Schechter approach to non-Gaussian density fields 




(51) 
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Using the spherical isothermal collapse model (Gunn & Gott 1972), it is possible to relate the 
mass M of a dark halo to the Lagrangian (pre-collapse) comoving length R (the sphere of radius 
R will give rise to an object that contains the mass M within the virialization radius) (e.g. White, 
Efstathiou & Frenk 1993), 

where R is in units of Mpc h _1 . Here V c is the circular velocity required for centrifugal support in 
the potential of the halo, given by M = /(WGH(z c )), where H denotes the Hubble constant, 
(e.g. Heavens & Jimenez 1999), and f c is the density contrast at virialization of the newly-collapsed 
object relative to the background. This is adequately approximated by f c = 178/Sl^ 6 (z c ) (e.g. Eke, 
Cole & Frenk 1996). 

In Figure 3 the dependence of (52) on the cosmology and the redshift of collapse is shown. 
It is clear that this dependence is very weak and can be ignored for our purposes. For a chosen 
power-spectrum shape and normalization, the dependence on the cosmological model of the PDF 
and the level excursion probability is therefore confined in S c (z c ). In other words, changing the 
cosmology is essentially equivalent to changing z c according to the right panel of Figure 3. 



Following the PS formalism, we can write 



TOO /*00 /'OO J\ 

P(>5 c \z c ,R)= d5 R P{5 R )= d5 R ' 

J Sc(Zc) J 6c(Zc) J —• 



-Xz c ) J8 c (z c ) J-oo 2vr 



e -iX8 R+ W(X) j (53) 



which, exchanging the order of integrations and integrating over S R , yields the exact and general 
expression 

i poo J\ 1 

P(>S c \z c ,R) = —J — exp[-iX5 c (z c )+W(X)} + -. (54) 

This is the exact expression for the probability that at redshift z c , the fluctuation on scale R 
exceeds some critical value <5 C . From this equation it is possible to obtain an approximate expression 
for the cumulative probability P(> S c \z c ,R) by expanding W(X) of eq. (27) to a given order in e. 
On the other hand, from knowledge of the cumulants of the cosmological density field smoothed on 
scale R, up to some order, it is possible to obtain an approximate expression for the level excursion 
probability P(> S\R). 

The integral (54) can be solved analytically by using the saddle-point technique as in (Fry 
1986) and (Lucchin & Matarrese 1988). We first perform a Wick rotation, A — ► iX, in the complex 
A plane, to get rid of the oscillatory behaviour of the integrand. This gives 

i pioo j\ i 

P(> S c \z c , R) = ^~ -r exp l-X5 c (z c ) + M(X)] + - , (55) 
2m J_ ioD X 2 



with M(X) = W(— iX). Next, let us introduce the effective action G(S e g), defined as the Legendre 
transform of M, namely G(S e s) = X5 c g — M(X), with S e g = dM/dX. From the definition of M 
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one has 

P(> 5 c \z c ,R) 



1 

2ri 



G'=ioo 
G'=-ioo 



d^cs^j^l exp [(5 cS - 5 c {z c )) G'(5 cS ) - G(S eS )] + \ 



G'(5, 



cS) 



(56) 



For large thresholds, 5 c (z c ) » 1, the above integral is dominated by stationary points of the 
exponential. These occur at G"{5 c r )(5 c (z c ) — 5 c d) = 0, i.e. at 8 e d = 5 c (z c ), since G"(5 c f[) = 
{(fM/dX 2 )' 1 > [as M must be a convex function of its argument (e.g. Fry 1985)]. A saddle- 
point evaluation of this integral then gives (Lucchin & Matarrese 1988) 



P{>5 c \z c ,R) 



1 (G") 1/2 
v^F G' 



exp(-G) 



S e B=S c (z c ) 



(57) 



Note that the remaining 1/2 term on the r.h.s. of equation (56) has been exactly canceled by the 
pole residual of the integral at G'(8 e g ) = 0. The latter formula provides an accurate approximation 
of the level-crossing probability P(> S c ), provided the cumulant generator M(5 e s) is analytic for 
all finite values of its argument, which is indeed the case for our approximated expression in eq. 
(50), but certainly not for the exact form of eq. (27). The validity of our last result will then rely 
on the smallness of the non-Gaussianity parameter e. It is important to notice here that since - as 
it is clear from Figure 1 - the skewness parameter is not small everywhere, we expect the PDF to 
be very sensitive to small e at least for model A. 



G(5, 



The approximate form of the effective action resulting from Legendre transforming eq. (50) is 

S3,R 



•cS) 



24 












T2 





off 



(58) 



The level crossing probability P(> 5 C ) to second order in e can be obtained from (57), with 
(57) and the following substitution: 



(G"') 1/2 
G> 



^eff 



2 

1 + ^ (72/3(0 2 )4 " « 2 /3<5 c 2 ff4^ 2) + 12a/4 2) + SS^ aS 2 eS a 2 R - 25 2 s a 2 R aSP) 



To linear order in e the effective action reads 



G(<5 ef r) 



24 



S 3 . 



R 



3 



(59) 



(60) 



Note that this function is only convex for <5 e fr < 1/^3 A> so * na * our resur t can be consistently 
applied only as long as 5 c (z c ) < 1/S^ R . 

We finally obtain 



P(> 5 c \z c ,R) 



OR 



2tt S c (z c ) 



exp 



l<fe) 



4 



, S 3,R X , \ 
1 5— Oc{Zc) 



(61) 
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Note that, for the Gaussian case, S^r = 0, this formula corresponds to the well-known asymptotic 
behaviour of the complementary error function valid where 5 C S> 2a. In particular, for S^r = 0, 
and z = this approximation introduces an error of 20% at R ~ 10 Mpc hr 1 and an error smaller 
than 5% for R > 20 h' 1 Mpc. 

For 5 c (z c ) > l/S^ij, i.e. for R < 10 Mpc and/or e > 10 -3 for model A, the integral (54) has 
to be performed numerically. To first order in e eq. (54) becomes: 

1 1 r°° d\ 

P(> S c \z c ,R) = --- — exp(-A 2 4/2) sin(A,5 c + A 3 /x 3 , R /6) . (62) 

2 7T J A 

In Figure 4 we show the result of the numerical integration for model A for masses in the range 
10 8 - 1O 12 M , redshift z c = 6, 8, 10 and e from top to bottom e = 10~ 2 , 5 x 10" 3 , 2 x 10~ 3 , 10~ 3 , 
5 x 10~ 4 , and the Gaussian case (e = 0, solid line). The main effect of the presence of a small 
non-Gaussianity such as e = 5 x 10~ 4 , is to amplify P(> 6 C \ z c , R) by a factor of order 10. 

For model B instead, as shown in Figure 2, the skewness parameter for the density field is 
small and has the opposite sign of e. For this reason the saddle-point approximation (60) works 
remarkably well on galaxy scales for |e| < 500, if z c < 10. Figure 5 shows the excursion set 
probability for model B for masses in the range 10 8 - 1O 12 M , redshift z c = 6,8,10 and e from 
top to bottom e = -300, -200, -100,0, 100,200. The thick solid line is relative to the Gaussian 
case (e = 0); curves above that have e < 0, below have e > 0. The fact that |e| ~ 100 is needed 
to create any sizeable departure from the Gaussian P(> 5 C \M), has an important consequence: 
all conventional inflationary models induce deviations from Gaussianity such that the skewness 
parameter for the peculiar gravitational potential 5^$ is bound to be \S^ ; ^\ < 10 (Gangui et al. 
1994); since ~ 6e, we conclude that, in the context of single-field inflation models the level- 
crossing probability P(> S C \M) is indistinguishable from the Gaussian one at least on clusters scale 
and below. It is important to stress that for this model, using values of e of order unity or even 
larger does not imply the breakdown of the perturbation expansion upon which our results rely, as 
the coefficients of that expansion, such as the skewness parameter, keep small even for |e| < 100. 

The most complex issue is that of finding a suitable expression for the function P(M\M') in 
equation (51) in the general non-Gaussian case and for generic filters. In this context, it is useful 
to define a 'fudge factor' / through the equation 

P{M\M') = j 0(M' - M) , (63) 

with 6 the Heaviside step function. As shown by Nagashima & Gouda (1997), in the Gaussian case, 
with sharp-k-space filtering, P(M\M') reduces to the conditional probability P(5m > 8c\$M' = dc), 
which is easily obtained from a bivariate Gaussian, using Bayes theorem; from this one immediately 
gets / = 2, as expected (e.g. Peacock &: Heavens 1990). In the non-Gaussian case (and/or for 
different types of filter) the value of / should be obtained from its very definition. Recently, 
Koyama, Soda & Taruya (1999) showed that in generic non-Gaussian models one can still write 

-1 



f(M, M') 



/ duj p(5 M = v \$m> = $c) 

.J S c 



(64) 
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with p(Sm = w\8m' = <5c) ~ p(5m = w — <5c), provided M' 3> M. Obtaining a similar relation for a// 
scales M' > M is a more complex issue, which would require a separate analysis 8 . Nonetheless, it is 
extremely reasonable to expect that this approximation works well at least in so far as the deviation 
from Gaussianity is weak. One then immediately gets the simple result 9 / ~ 1/P(> 0|M). In our 
case, one can compute the cumulative zero-crossing probability P(> 0,M) numerically, starting 
from eq. (62). 

The comoving mass-function of halos formed at redshift z c can thus be immediately obtained 
by differentiating the integral equation (51). This allows the standard PS formula to be extended 
to the non-Gaussian case in the simple form (Lucchin Sz Matarrese 1988; Colafrancesco et al. 1989; 
Chiu et al. 1998; Robinson et al. 1999a,b; Koyama et al. 1999): 



8irGM 



dM 



(65) 



where we implicitly assumed that the value of / has negligible dependence on M, The factor / should 
account for the physical constraint that all of the clustered matter in the Universe must be included 
in bound objects of some mass [i.e. J °° dMMn(M, z c ) = 3-ffgf2om/87rG] and simultaneously solve 
the cloud-in-cloud problem. 

For the non-Gaussian models considered here, the factor / is not much different from the 
Gaussian value, provided e is small, i.e. for small departures from the Gaussian behaviour. In 
Figure 8 the fudge factor / ~ 1/P(> 0\M) is plotted as a function of the skewness parameter, 
using a top-hat filter in real space. The useful range is on the left of the vertical dotted line, that 
is M > 2 x 1O 1O M for model A with e = 1CT 3 and M < 4 x 10 15 M Q for model B with e = -100. 
We can therefore conclude that, for mild non-Gaussianity, the correction to the usual fudge factor 
/ = 2 is negligible when compared to the effect of the skewness term in the exponential, justifying 
therefore the use of / = 2 throughout. 



4.1. Model A 

In model A the skewness parameter S^r is scale-independent on large scales, therefore, as a 
first step, we can neglect the dependence of S^r on mass. In this case, in the above differentiation, 
we arrive to a very simple generalization of the PS formula, which has the same functional form, 



8 In the Gaussian case, and for sharp-k-space filter, one finds p(5m = u\6 M ' = S c ) = p(<5jj = u>—6 c ), for all M' > M, 
where M is defined by a 2 (M) — a 2 (M) — a 2 (M'). Inserting this result in eq. (64) yields / = 2, independently of the 
mass (Nagashima & Gouda f997). 

9 Note that, if the non-Gaussian PDF depends on the filtering mass only through the variance, i.e. P(Sm) = 
a~^^(S M /cFM), as assumed by Willick (1999), then P(> 0|M) = P(> S C \M = 0), which yields an alternative 
expression for /, used by some authors (e.g. Robinson & Baker 1999). In the most general case, however, the latter 
result is not valid. 
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provided one makes the replacement 

5 c (z c ) -> 5 c (z c ) 



S3,R 



1/2 



= 6*(z c ) 



and therefore: 



n(M, z c ) ~ / 



3-Hq ilom 



<5*(^c) 



8vrGM 2 



■ exp 



Sj(z c ) 



din ou 



d\nM 



(66) 



(67) 
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This simple result is completely analogous to that found by Lucchin & Matarrese (1988), for 
hierarchical- type statistics: the effect of the deviation from primordial Gaussianity is that of shifting 
the PS characteristic mass to higher or lower masses, depending on whether the skewness is positive 
or negative. Note that the effect of the non-Gaussian tail becomes more and more important at 
higher redshifts. As noticed by Fry (1986), 'the departures from Gaussian can be appreciable while 
fluctuations are still small'. The assumption that S3 .r is scale- independent is correct on relatively 
large scales, where R > 20 h^ 1 Mpc, corresponding to masses M > 1O 14 M . Substitution (66) 
applied to the standard PS formalism agrees with the numerical results on scales R > 20 h^ 1 
Mpc, for e < 10~ 2 . The dependence on the cosmological model is confined to S c (z c ), as discussed 
before, once the present-day power-spectrum shape and normalization have been chosen, and the 
normalization i/gl^om- 

However from Figure 1 it is clear that the scale dependence of S^ R is not negligible for galaxy 
masses, moreover for R < 10 i.e. M < 10 13 the saddle-point approximation does not hold. Numer- 
ical evaluation of (54) yields dP(> 5 c )/dM; the corresponding n(M,z c ) is shown in Figure 6 for 
the same range of masses, redshift of collapse z c and non-Gaussianity parameter e as in Figure 4, 
assuming / = 2, according to the discussion above. 



4.2. Model B 



It is clear from the scale dependence of the skewness parameter that, for departures from 
Gaussianity of the kind of model B, clusters scales are better than galaxy scales to probe primordial 
non-Gaussianity. It is worth to notice here that also for model B in the absence of the transfer 
function, S^^r is scale independent and indeed the skewness parameter becomes scale-independent 
on very large scales (R > 100 h^ 1 Mpc), where the transfer function is unity. On clusters scales 
Ss t R is not scale-independent, but for e < 200 and redshift z c < 2 the saddle-point approximation 
is still valid. For the comoving mass-function of halos formed at redshift z c we obtain: 



n(M, z c ) ~ / 



3-Hq f^Qm 



1 



8ttGM 2 v^Fctm 



8 2 c (z c ) 



dSs )M . 5*{z c ) d\na M 



+ 



2 3^/1- S 3 ,mS c (z c )/3 din M a M dlnM 



exp 



(68) 

Also in this case the dependence on the cosmological model is confined to 5 c {z c ) and in the nor- 
malization H^om. Figure 7 illustrates the effect on the mass- function on clusters scales for the 
non-Gaussianity of model B. For two different redshifts of collapse (z c = 1 and 2) the ratio of 
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the mass-function for model B to the mass-function for a Gaussian field [n(M, z c ) / riG SLVL (M , z c )\ is 
plotted as a function of the mass in solar masses. Also in this case we can assume / = 2. The choice 
for the e parameter is, from top to bottom, e = —100, —50, —10. It is clear that for high masses 
one is probing the tail of the distribution, that is most sensitive to departures from Gaussianity. 
The forthcoming X-ray Multi-Mirror (XMM) galaxy cluster survey (Romer et al. 1999), will allow 
the number density of clusters to be accurately measured even at z > 1. In particular the XMM 
cluster survey will provide a complete survey of clusters with masses M > 1O 14 M at z < 1.4 over 
800 square degrees. Preliminary calculations suggest that it will be therefore possible to detect 



5. A worked example: galaxies at z > 5 

The technique presented in the previous section allows one to measure accurately the amount 
of non-Gaussianity on a given scale. Although the traditional route has been to use the abundance 
of clusters, we will illustrate our technique by using the observed abundance of high-z galaxies 
(e.g. Cavaliere & Szalay 1986; Kashlinsky k, Jimenez 1997; Peacock et al. 1998). There are some 
advantages in using high-z galaxies: they sample directly the galaxy scale and the objects are always 
inside virialised halos. One major disadvantage is that the mass is not accurately determined. On 
the other hand, as larger telescopes, such as NGST, get on line, it will be possible to obtain volume 
limited samples of of high-z galaxies and from high-resolution spectra, dynamical masses will be 
determined more accurately. As seen from eq. (66), it is clear that the higher the redshift the more 
sensitive the PDF is to non-Gaussianity. Up to date there are 6 galaxies with spectroscopic redshifts 
observed between z > 5 to z ~ 7 (Dey et al. 1998; Spinrad et al. 1998; Weymann et al. 1998; Hu 
et al. 1999; Chen et al. 1999) and in Table 1 we have compiled their most relevant characteristics. 
These high-z galaxies have been found in relatively small areas of the sky: for the purposes of this 
worked example we will use a scanned area of 3 x 10~ 4 square degrees for each galaxy 10 . 

As pointed out before, the masses of this high-z sample are unknown, but their current obscured 
star formation rate is known. The first step is to correct the observed star formation rate for the 
presence of dust. In (Jimenez et al. 1999) a case has been made for the correction factor for the star 
formation due to the presence of dust in a starburst being ~ 5 (independently of the magnitude 
of the host galaxy). This is confirmed by the independent arguments presented in Peacock et al. 
(1999). The corrected star formation rate is presented in Table 1. Since the age of the Universe 
at z ~ 6 is already about 0.7 Gyr even for EdS, assuming the present age to be 13 Gyr, it is 
very unlikely that galaxies have been observed just when their stellar populations are born. Note 



1 'Four of the six galaxies considered here have been found in the Hubble Deep field, where the total area scanned 
is about ~ 1/800 of a square degree. The other two galaxies have been found in Keck observations, where the area 
scanned for each of them is about 7.4 x 10~ 5 square degrees. For this worked example we will therefore make the 
conservative choice that there is on average one galaxy per 3 x 10 -4 square degrees. 
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that galaxies are likely to be dust-enshrouded for about 0.02 Gyr (Jimenez et al. 1999) and thus 
be invisible for about 3% of the Hubble time at z ~ 6. The star formation rate of the galaxy in 
the past is more likely to be higher (Heavens & Jimenez 1999) than the present observed one. To 
account for this we will then assume that the galaxy has been observed when it is 0.1 Gyr old 
(i.e. about 10% of the Hubble time at z « 6) and has been forming stars at the current (constant) 
observed rate. Using this, we can estimate what the mass in stars in the galaxies (see Table 1) and 
in dark matter is, using a simple isothermal profile for the halo, i.e. M = /(10GH(z c )), where z c 
is the redshift of collapse 11 and assuming that Oob = 0.015/i~ 2 . Following Peacock et al. (1998) we 
will carry out our analysis in terms of a v . For an isothermal halo a v = V c /y/2. Table 1 shows that 
a v = 70 — 200 km s _1 . The number density of high-z galaxies is N(Q = 1) > 3.6x 10 _3 (/i _1 Mpc) -3 
[N(Q = 0.3, A = 0.7)> 8.3 x lO^/i^MpcT 3 ]. 

To compare these observations with our theoretical results we integrate the mass-function to 
obtain the number of galaxies observed per unit volume with redshift of formation 6 < z c < 8 
and masses 2 x 1O 1O M < M < 4 x lO n M . It is worth noting that our mass estimates slightly 
underpredict the masses that are derived using the width of the Lya line (Dey et al. 1998). Gaussian 
initial conditions give N(fio = 0.3, Ao = 0.7)= 5.2x 10 _5 (/i^ 1 Mpc) -3 , which is about a factor sixteen 
lower than the observed value. For Model A with e = 10 -3 , we obtain N(S7o = 0.3, Ao = 0.7) = 
1.3 x 10~ 3 (/i _1 Mpc) -3 , which is much closer to the observed value. As we already noted before, a 
very small non-Gaussianity parameter (i.e. a very small departure from the Gaussian behaviour) 
has a dramatic effect on the statistics of high redshift objects. This is, of course, very preliminary. 
For the tail of the PDF, the mass function is very sensitive to the overdensity threshold, the mass 
determination and the redshift of formation. In this example we have assumed that the redshift of 
formation is the redshift at which the object has been observed. This is of course quite a conservative 
assumption, yielding a lower limit on the non-Gaussianity parameter. It has been shown that the PS 
algorithm with a fixed threshold 5 C does not reproduce very accurately N-body simulation results 
(e.g. Sheth & Tormen 1999). However, the deviations on the mass scales considered here are 
never larger than a factor of a few, well below the uncertainty due to the mass determination (see 
below). PS predictions and N-body results can be reconciled by making the threshold for collapse 
scale-dependent (e.g. Sheth & Tormen 1999; Sheth, Mo & Tormen 1999). This will change the 
absolute prediction on the number density of objects but will not strongly affect the ratio between 
the Gaussian and non-Gaussian number density predictions. The spherical collapse assumption 
underlying the PS algorithm is the main reason why theory does not reproduce accurately N-body 
results: the introduction of the ellipsoidal collapse improve sensibly the agreement (Lee & Shandarin 
1998). In a subsequent paper we will address the issue of triaxial collapse for non-Gaussian initial 
conditions. 



11 Although we have assumed that the stellar population of the galaxy has to be about 10% of the Hubble time at 
the observed redshift, we will conservatively consider that the dark halo has just collapsed at the observed redshift. 
This is what is expected in hierarchical models, where the stellar population is generally made in previous generations 
inside smaller haloes. 
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Maybe the most pressing issue is that of the mass determination of high-redshift objects. For 
example, if the mass estimate for the galaxies in the worked example is wrong by a factor of 4, the 
number density would change by a factor of 20 (at z c = 6), and agreement between predictions and 
observations could be obtained without resorting to e ^ 0. 

Indeed, better observations will improve the determination of the mass of the dark halos and 
of the number density of high redshift galaxies. A measurement of e and thus a determination 
of the amount of primordial non-Gaussianity (if any!) present on galaxy scales will then become 
possible. 

6. Conclusions 

We have shown that it is possible to test whether initial conditions were Gaussian, by looking 
at the abundance of galaxies and clusters at high redshift. Small deviation from Gaussianity have 
a deep impact on those statistics which probe the tail of the distribution: for this reason the 
abundance of high redshift objects such as galaxies at z > 5 and clusters at z > 1, is very sensitive 
to small departures from the Gaussian behaviour. 

We used a parameterization of primordial non-Gaussianity that covers a wide range of physi- 
cally motivated models: the non-Gaussian field is given by a Gaussian field plus a term proportional 
to the square of a Gaussian field, the proportionality constant being the non-Gaussianity param- 
eter. We applied this model to the primordial density fluctuation field and to the gravitational 
potential fluctuation field, and we assumed that the departures from the Gaussian behaviour are 
small. Given this parameterization of non-Gaussianity, we were able to obtain an analytic expres- 
sion for the probability density function (PDF) of the smoothed density field. This approach is 
somewhat different from what is found in the literature (e.g. Robinson et al. 1999a, b), since we do 
not assume a particular form for the non-Gaussian PDF, but we derive it, given a parameterization 
of the primordial non-Gaussianity. 

We then introduced a generalized version of the PS approach valid in the context of non- 
Gaussian initial conditions, also tackling the cloud-in-cloud problem. This enabled us to relate, 
analytically, the non-Gaussianity parameter to the number of high-redshift objects. The strength 
of our method is that we are able to properly take into account the smoothing scale, or mass, 
dependence of the probability distribution function, and that our results are analytic. In looking 
at the likelihood of rare events for a non-Gaussian density field, one is probing the tail of the 
distribution, and analytical results are not only elegant but also extremely important. 

The main technical results of the present paper can be summarized as follows: a) Eq. (27) 
is the exact analytic form of the cumulant generator for our non-Gaussian models: it can be used 
to generate cumulants of any order by simple differentiation, b) Eq. (54): this equation gives an 
exact and general analytic expression for the cumulative probability that, at a given redshift z c , the 
fluctuations on scale R exceed some critical value 5 C . From this expression it is possible to obtain 
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the cumulative probability to any order by expanding the cumulant generator (36) with (40) and 
(41) as we did in (50). c) Eq. (61): it is an analytic expression for the cumulative probability valid 
for small deviations from the Gaussian behaviour, d) Eq. (68), with the definition (66): it is an 
analytic expression for the comoving mass-function of halos obtained within the PS approach. For 
a given power spectrum the dependence on the cosmological model is straightforward: it is confined 
to the normalization factor H^Clom and the threshold <5 C . 

Still an important issue remains unsolved: how to accurately determine the mass and the 
formation redshift of observed high-redshift objects in order to put tight constraints on primordial 
non-Gaussianity. This issue will be addressed in a forthcoming paper. 

LV acknowledges the support of a TMR grant. LV and RJ thank A. Berera, A. Heavens and 
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Figure 1: Model A variance, skewness (/ig R ) and skewness parameter (S 3 R ) as a function of the 
radius of the top-hat window in real space and of the mass. To produce this plot we assumed: a = 1, 
a scale-invariant initial power-spectrum P^ oc k, CDM transfer function, f^om = 0.3, normalized to 
produce a%(z = 0) = 1 (more details in the text). 

Figure 2: Model B variance, skewness (/^r) an d skewness parameter (S^ R ) as a function of 
the radius of the top-hat window in real space and of the mass. To produce this plot we made the 
same assumptions as for Figure 1 except that now P§ oc A; -3 . 

Figure 3: Left panel: The redshift dependence of the threshold for collapse A c in different cos- 
mologies. It is clear that the dependence on redshift and on cosmology is weak. This justifies to 
assume A c =1.687 throughout. Middle panel: Dependence of the relation between dark mass and 
radius on the cosmology and the redshift of collapse z c . The solid lines denote an Einstein de Sitter 
Universe, the dashed lines a flat model with fiom = 0.3 and Ao = 0.7 and the dotted lines an open 
model with flo m = 0.3. The dependence is very weak and can be safely ignored. Right panel: The 
redshift dependence of the inverse of the linear fluctuation growth for different cosmologies. 

Figure 4: Model A P(> 6 C | z c ,R) for z c = 6 (left panel), z c = 8 (center), z c = 10 (right panel) 
and e = 10~ 2 , 5 x 10~ 3 , 2 x 1CT 3 , 10 -3 , 5 x 1CT 4 . The thick solid line is relative to the Gaussian 
case. It is clear that the main effect of the presence of a small nonzero e is to boost P(> S c \ z c , R) 
by at least a factor ~ 10. 

Figure 5: Model B P(> 6 C \ z c ,R) for z c = 6 (left panel), z c = 8 (center), z c = 10 (right panel) 
and e = -300, -200, -100, 0, 100, 300 (from top to bottom). The thick solid line is relative to the 
Gaussian case. It is clear that | eb | of model B need to be S> of model A to show a noticeable 
difference from the Gaussian case. Moreover es < is needed in order to enhance the non-Gaussian 
tail of P(> 5 C | z c ,R). 

Figure 6: Model A: The comoving mass-function of halos formed at redshift 6 (left panel) 8 
(center), 10 (right panel), for different values of the non-Gaussian parameter e, as in Figure 4; from 
top to bottom: e = 10" 2 , 5 x 10~ 3 , 2 x 10~ 3 , 10~ 3 , 5 x 10~ 4 . The thick solid line is relative to the 
Gaussian case. 

Figure 7: Model B: The non-Gaussian effect on the mass-function on clusters scales. For two 
different redshift of collapse (z c = 1 left and z c = 2 right) the ratio of the mass-function for model 
B to the mass-function for a Gaussian field [n(M, z c )/nQ au (M, z c )\ is plotted as a function of mass. 
The choice for the e parameter is, from top to bottom, e = —100,-50,-10. It is clear that for 
high masses one is probing the tail of the distribution, that is most sensitive to departures from 
Gaussianity. 

Figure 8: The fudge factor / ~ \/P{> 0|M) as a function of the skewness parameter. The value 
/ = 2 is exact for the Gaussian case with sharp-k-space filter; the solid line shows the value of / 
as a function of the skewness parameter S^r. The useful range is on the left of the vertical dotted 
line, i.e. M > 2 x 1O 1O M for model A with e = 10~ 3 and M < 4 x 1O 15 M for model B with 
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e = —100. We can conclude that, for mild non-Gaussianity, the correction to the usual fudge factor 
is negligible, justifying therefore the use of / = 2 throughout. 
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Table 1: Main properties of observed galaxies at z > 5. Column 1: [1] Dey et al (1998); [2] Spinrad 
et al. (1998); [3] Weymann et al. (1998); [4] Hu, McMahon & Cowie (1999); [5] Chen, Lanzetta 
& Pascarelle (1999). Column 2: spectroscopic redshift of the galaxy. Column 3: area of the sky 
(in square degrees) surveyed per galaxy. Comumn 4: observed star formation rate in solar masses 
per year (assuming h = 0.65). Column 5: dust corrected star formation rate (see text). Columns 
6, 7: estimated mass in stars and dark matter (see text). Column 8: velocity dispersion obtained 
assuming an isothermal halo (see text). 



Ref. 


z 


deg 2 


obs.SFR 


corr.SFR 


est.M* 


cst.M dark 


<j v 




(Moyr-i) 


Mq 


Mq 


(kms -1 ) 


[1] 


5.34 


7.4 x 10~ 5 


4 


20 


2 xlO 9 


2 xlO 10 


70 


[2] 


5.34 


3 x 10" 4 


13 


65 


6.5 xlO 9 


6.5 xlO 10 


100 


[2] 


5.34 


3 x 10~ 4 


13 


65 


6.5 xlO 9 


6.5 xlO 10 


100 


[3] 


5.60 


3 x 10~ 4 


8 


40 


4xl0 9 


4 xlO 10 


86 


[4] 


5.74 


7.4 x 10~ 5 


40 


200 


4xl0 10 


4xlO n 


187 


[5] 


6.68 


3 x 10~ 4 


40 


200 


4xl0 10 


4 xlO 11 


200 
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Fig. 2.- 
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Fig. 6.— 
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Fig. 7.- 
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